Parameters

Document rendering

Analysis input/output

input_folder <- "raw_input" # Where all the large input files are. Ignored by git. 
output_folder <- "results" # Where plots will be saved
output_format <- "pdf" # The file format of saved plots
pub_fig_folder <- "publication"
revision_n <- 1
result_path <- function(name) {
  file.path(output_folder, paste0(name, ".", output_format))
}
save_publication_fig <- function(name, figure_number) {
  file.path(result_path(name), paste0("revision_", revision_n), paste0("figure_", figure_number, "--", name, ".", output_format))
}

Analysis settings

min_p_value <- 0.05

Introduction

Dependencies

This analysis requires Bioconductor packages to work. These packages are not available on CRAN and must be installed using the following code. The code below is not run during the compilation of this page, but is provided so others can quickly find the right software to replicate this analysis.

source("https://bioconductor.org/biocLite.R")
source("http://bioconductor.org/workflows.R")
workflowInstall("rnaseqGene")
biocLite("GO.db")
biocLite("org.Hs.eg.db")
biocLite("airway")
biocLite("DESeq2")

Read in data

Starting from SummarizedExperiment

http://www.bioconductor.org/help/workflows/rnaseqGene/

library(rnaseqGene)
library("airway")
data("airway")
se <- airway
se$dex <- relevel(se$dex, "untrt")
se$dex
## [1] untrt trt   untrt trt   untrt trt   untrt trt  
## Levels: untrt trt
library("DESeq2")
dds <- DESeqDataSet(se, design = ~ cell + dex)

Pre-filtering the dataset

nrow(dds)
## [1] 64102
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)
## [1] 29391

The rlog transformation

rld <- rlog(dds, blind=FALSE)
head(assay(rld), 3)
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003   9.385683   9.052608   9.516875   9.285338   9.839085
## ENSG00000000419   8.869616   9.138271   9.036116   9.075295   8.972126
## ENSG00000000457   7.961369   7.881385   7.824079   7.921490   7.751699
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003   9.530311   9.663255   9.277699
## ENSG00000000419   9.131824   8.861534   9.060905
## ENSG00000000457   7.886441   7.957121   7.912123

Differential expression analysis

dds <- estimateSizeFactors(dds)
dds <- DESeq(dds)
(res <- results(dds))
## log2 fold change (MAP): dex trt vs untrt 
## Wald test p-value: dex trt vs untrt 
## DataFrame with 29391 rows and 6 columns
##                    baseMean log2FoldChange      lfcSE       stat      pvalue
##                   <numeric>      <numeric>  <numeric>  <numeric>   <numeric>
## ENSG00000000003 708.6021697    -0.37415193 0.09884432 -3.7852648 0.000153545
## ENSG00000000419 520.2979006     0.20206144 0.10974240  1.8412340 0.065587276
## ENSG00000000457 237.1630368     0.03616620 0.13834538  0.2614196 0.793768939
## ENSG00000000460  57.9326331    -0.08446385 0.24990676 -0.3379815 0.735377161
## ENSG00000000938   0.3180984    -0.08413904 0.15133427 -0.5559814 0.578223585
## ...                     ...            ...        ...        ...         ...
## ENSG00000273485   1.2864477     0.03398815  0.2932360  0.1159071   0.9077261
## ENSG00000273486  15.4525365    -0.09560732  0.3410333 -0.2803460   0.7792120
## ENSG00000273487   8.1632350     0.55007412  0.3725061  1.4766847   0.1397602
## ENSG00000273488   8.5844790     0.10515293  0.3683834  0.2854442   0.7753038
## ENSG00000273489   0.2758994     0.06947900  0.1512520  0.4593591   0.6459763
##                       padj
##                  <numeric>
## ENSG00000000003 0.00128686
## ENSG00000000419 0.19676183
## ENSG00000000457 0.91372953
## ENSG00000000460 0.88385059
## ENSG00000000938         NA
## ...                    ...
## ENSG00000273485         NA
## ENSG00000273486  0.9062268
## ENSG00000273487  0.3389275
## ENSG00000273488  0.9039857
## ENSG00000273489         NA

Converting ENSEMBL IDs to GO IDs

library(GO.db)
library(org.Hs.eg.db)
res$go_id <- mapIds(org.Hs.eg.db,
                    keys=rownames(res),
                    column="GO",
                    keytype="ENSEMBL",
                    multiVals="first")

Fitlering results

Remove genes with no GO annotation

nrow(res)
## [1] 29391
res <- res[!is.na(res$go_id), ]
res <- res[res$go_id %in% keys(org.Hs.eg.db, keytype = "GO"), ]
nrow(res)
## [1] 15351

Remove insignificant genes

nrow(res)
res <- res[(! is.na(res$padj)) & res$padj <= 0.05, ]
nrow(res)

Remove genes with small changes

nrow(res)
res <-  res[abs(res$log2FoldChange) >= 0.5, ]
nrow(res)

Getting classificaiton

term_class <- function(x, current = x, all_paths = TRUE, type = GOCCPARENTS, verbose = TRUE, 
                       valid_relationships = c("is_a")) {
  # Get immediate children of current taxon
  parents = tryCatch({
    possible_parents <- as.list(type[x[1]])[[1]]
    if (! is.null(valid_relationships)) {
      possible_parents <- possible_parents[names(possible_parents) %in% valid_relationships]
    }
    names(AnnotationDbi::Term(possible_parents))
  }, error = function(e) {
    c()
  })
  
  # only go down one path if desired
  if (! all_paths) {
    parents <- parents[1]
  }
  parents <- parents[parents != "all"]
  
  if (is.null(parents)) {
    return(c())
  } else if (length(parents) == 0) {
    cat(length(x))
    return(paste0(collapse = "|", AnnotationDbi::Term(x)))
  } else {
    next_x <- lapply(parents, function(y) c(y, x))
    
    # Run this function on them to get their output
    child_output <- lapply(next_x, term_class, all_paths = all_paths, type = type)
    output <- unlist(child_output, recursive = FALSE)
    
    return(output)
  }
}
cc_class <- lapply(res$go_id, term_class, all_paths = FALSE, type = GOCCPARENTS)
## 55344646434566443568574246946246485566422686555546464854425646355675264568667442845484545445734645543554335646524765144556452444453655326645667643665636437266365553666653554443547755544665656656617564665446764556665646573766647856634425556465333555454566443555466577663445647654765543454444555234434445463554664643554674427516626436466345327355425644457463455625334565764765555635665654446226564456762536437552426424656415646765565365355455756667445245644673555633443446654543448365557356745464253446355646785353564652445438566974644544762464564634664566546216445464625664546666543764666444444557654654656476665655364555668566664854733555844645856664356564664764543466574444544356564547444262664234462626633544766644566543556647578368774345366535665564555456356267443576546528764534646666576656356635466665625643444667455645576344444653635446345333645354667767555366575455346454433732665346654666546466456345545645366544686434676466536673255474444236441454645576557246444556354557336765458552565333524444754444253455444543665426466175443537538636645585636574355342552445655557566345555345644757445546655644434556167755254476552563345555354668456653657854564542555645545433556535564335646574743483668645444465355767753557465646345643455555652556456665565755459696647646565456448454854545446325563564144555435463474564654668563646663464436664536258563553454463443634557553564532442563665347464855535674557465267425626544442455633631664266666646655566654445644646592354355636446452525667273554452465531665644744355765545245442346544334444665554464753645566524645566676564644467265355623634565665426463954554635374633443685444655654624566442445442636345338566874255665749365544663266446546666434376545635875447568575244665554476457534564445646656432332565556423465455545643363653654624424644655543523554424247634445736554565644665565434525245648462424664425174667444543645645663462442426354655334565444664664345656544445442426456554333534235555445226557466427453466452355885266552564543866754823454536276665554535554328424454455474453445344555544572383854525444254153554469572265563545646326556565635354526685344274756384355575133565246443443535225577464654626645465944665455364555543266254375454664425663635469256464555532664444645645153766643556454444455213625466563353656664635452446463434537564635565452645654333346653543455445657646648341453536557355544554654362364544653566246466565334356655246424446454654264466554456555545366546144356265546659564564554664365453424456453255736464354554425234534455553474810335445455444354445764674435344375465246644535444456466236654536635334664464443454454232666452343456652645625353467445555754233555632495655145235443375654634425465473554665474555567446467444434572269427544343443444445745544442554547452353452442454654546435454522444656777254345444484254243465456684585425564554564535555545455253522445646654547865735452466566494446423445535954153345756533463454334666446344623455424455341464366457726846646264654456624635344656744675445443455355534656456354544344456532443664576614662544234526654544352447545553454555342646345656361665525656534636555546456442545444434484446554425464444625424427552442324743546245485528454425234255664255444345462456544567645347654444573556324634532434474363444747733656444824445575645755622552637232764854445543453562454242154
mf_class <- lapply(res$go_id, term_class, all_paths = FALSE, type = GOMFPARENTS)
## 2736379397755347564885483563353810314333534101036933532757107655515656196678179553455617265757354536333537735646714155459415875346738338337314723336785511874353433359534315935382751553368656333453892334510878877610465527581535623638566535736533756366753754314166363888553525779810435105352955235533466433353383313333363938335747731579955636683955365735474768753466576986681542347365857633943243597653115335737635535483566335636565337355165610941431033832916349485766339951673896476106366910591136875410656342749665645576324646873286656264196585693654837563478106575384476589574336641067334854536951073616431553310134588384486652766543566455483673545564657374433868255533575476867653653748593188335510466363334236851038366665441016653848347684849555941065547866558385335373363337107165367738636346618654566744764875366538431083910445753633533385636310737759455386158433995510363865356103374579645564105365862553618364835234134433446851466634684633388551593763656941463364554655744369546447545355106337346563763365635672658433631543710633310551453513134536832853637657643918799445333335368666521432385434255858610144733763443368683368633335114653343863553710118173384936363454334464961145657568741645376156523923659610936113557633377455163456112745366735337634333235239110535356378376377413667254362194974373663379610471663355567234675545626333510346173756357916333337536556331034533753323747544457104366335599365635333733531695754571353344333395357331169435438966569561549531010633751510103675768765436655335103546385514753553164343857434941210111031473567386654681331713656710433533467910334533545833531068664753103564365710371035145455353638861086663653738833145525872346367556546313689315536533733174532679685985346713351453663876543133753453723646538333961513556567543556513610326977710456333235659356777687107534555843338176167666344346333186334333186435645635184355615183107242155810335863333333761436846785733338653363666676456633846236768661063363105102454546331062333638757596665266463377356756536723765859418779354333463531633644389145837535333553295453610103388168565336337763108534883374346491064666363534353635749855559553845395443854959166345366455565105563573633728835641443683333335674666543566523616375697171055561056935431589748655610432353359647567634533335334137310356336577833445510145855548333393531577834311331063339375116666133810313448934103733755453437635441538957513349149564394631375556410195686757767167935536552435331755643543105833455567345337797556683873639233714559111058677338577836186725104343554564661016674461354545755333567331565753535851378555531363336644399933351510933775310383474105543433654667157635443451037496558357373843433436556767555103387561464313343575363710851356556335556565856364566636869525448516338833417336363393375634688573315656577596634535833527365431031138551061334613695863357546446884333337353973973565335354345573765386833365553315355455186105635836543353574456635155355103353736669663375767355829256225354363851694553865731051746633333789578565181264535445635363385351153666338357438343337334343635997753367333565333563103546378583627711473635373331131773738584738535583372453131028466442753348338468573336635127762333877743734316276653765356333355344335910122563537358234883338776431737538108273645333561757776373612165137933357537555133644766109577571663461353141077733735105310333554451010451335167637145832654364549376434338963313644393558543461571037115674336481108737347488159957666688753731234636735576745691014676563435376353366355674648653887446661376433569167364175344263356338331111735381736654133343453134385343757686163697386587361437651333355136737776355145368356115471457331096655433843561556331551576710651061035333445558573387496279344353773797575337736851536548395511756617765656536333391686653813792311315555998856436355395373666645567534515335353336863462843310453392637453336386566938685634611110335338105218261066337413254338446858733346575353438381694576935166275365956839674656337246661654610533375376553756316453337535455676937337259334536343451336103835733126101254353544356359463854331053333653233663336315106351336718275335315334336393533531033345856753735363845345746155346534433138673810676691473355831171533563644356103433351765736358336731073599671573961364353677333559853836576510334468464474395131743395356533354253459598887497357396668463568654610237557453665663343535153653693337310563545649352631834365357777743323667103235243355343394378731977361423769328363533934383436855434363353143338833533466581485153561096636641594811066531065555133343567233361853564310543528223365739757931101043863561538463313556253786314368733353557173633313783333723338346554455312357341073713210243973645667526513174696535353564353437355555555331735177538354665813433164255719531037556635637741074348634343363445339385352144635374343536364566393434475663833353133117336436569855811745104742363438267313655635333693710331033433664777315653536107364756333566543733375151105433669331654577353643376510313356633937777667533337787855535645345444638653356333157553396637105545865752376553967645277537783433392445355178334773635343410488432335574538533562365233591152876933633683355786896563148333313644535355196573286756733766933513611634346415173164333977710453763351352331081551971145577483534653435336533137375533663616776684152338563510779833378754578346164631041075481764374383364486108733445573578496754557555584567493385553167441035555674154357935687752555553363571173632351363861637367310875634551033635535493258310839555733554815645743538523351075493775335335655635567561353455967731357355354433136595835855336444633674665484933373574555731016158386476356633519564103534555285539366544673973953357734353453510833755326562784552383101537353610333436546553610499339253154738545557845333731061355161011263373558661033115935853186104936868465399743595524593316179755345493346576173534565575593146765376855134658368754538155311027554567155373774733783553476355396387333759756955686361510567337348736865336959333676439115635841483359731635573736376553176536555434358173436575786456533655383105845649453535159973755394369375356381565547638310427867557273356339710331083163316337751061385553513119771633467553544973931375337845159437463513544656681103853557726551979133733363535343533595275646546535735843366541532555365745733363397104535215710371151035865733376551065783334521673386461397743351313773677734105633543510354533597357757345567193664465356661533534751665552335457355539953475313935911636556947555333536463855535313566537710553555233554555103355533109553663105105310158463333415553973311057571937583573131035106372471091105321746623538439535631386535795636434359587368733333155551353356655733485431341036385583983353993561310586466379337765555565543410943516974543128857343864491663391148346538636931036753554531833351105385313333637103132656103335735210771477191515468453573555593575836585235565775596375355426353953545435639445104733946671076653583653533348246316653376443343534333763517553565555665735535363578751531013295375555663577333538655334386564536566143353910477544536661374315653455756587341663552561055559595455155553102553513555758537557635615510855967365344385567653510310344555964653551037357233555553436131055105365555435735109107325673445553147137105537554176564510659597173533455357415953610533455557353585473589510555555235437463105810557357754105535561347111063818533553726745573359751955531085855535410617543673352755463521082574964543553775103537748377333553363335137363555658535446416753655437106515836184569445397454557617343349103193333551533447334552375335566666666553738535343964233373231815831736343355914135335637756746533355355114338531731555751878451517561336375105624588342453101110754535365753636132634107518433537373373975556576673367357410311034743455615153518334547554465387675865156644837373710175414576341557555510661317673525452104574356354365633371655534652541686535522673471356531486553356107934373413566339433375152456746565355165557610313513557654266310675667866466662566633425217333676757735351035537108335645355714684617311651566658338655535347495778636155
bp_class <- lapply(res$go_id, term_class, all_paths = FALSE, type = GOBPPARENTS)
## 8849745556574541097587561045114979797697776105748458776966744910457975976969101011999779856710966569741012988464127667889679657956669125119117411754679896554644410655968544511912866557971098951056789477686568741195696644851049475699954968875912895897541071212911510796757107485991159397712751141079994755791298481177109755571294879696666961161196910551094576412410568109101043119465575461049499111197977599677654787498791012597986944755885651010596771077411749839697410779951051191095610446469107569810574498779986797559799557101155555885859584295989554889107751076676939687741177851079471151111659751095910987810125910128410129879511557611959598758465574748761111511887849996465577554671057891074586109745499676781112671257911111041142711115559710764710981112655810510107511108757117479891169697957799589777991195799577996910115976879554841074866871167877961147767699498579712655551045585294971291191047555971167611545442599479549119577675955749854512966546105695115724127798994958899971079107114991194111157748911510654784891010756105101195581010791165912894595410125851164106779997799949795127129965891210104744961055876645410974554979111171298577106106106771059887111175894116134711976510496959121159781259610799126109115699441056449654610971158571110861241279108641157118410947461010779412119359559969696997556997117965775896564610576669411756991159751253751045648109101165911868810751011999106754746879965101255995125886558512996910105685924799766574581074884946558469119747949107769994795571056812997910669510511663984867664610856698129499555596875111176759977659105641067573959986666961111115711510794639109659912591096996128591291010645654669279610104951091165105910106989799459115774498595878456927610781277981059971061396610447711594115105954479664685695911956912810912118711949128119859446811995487109109779985859961191061398697118969864776697676757433795967669911796929510254777551095115597674676794774498104867114597991199114479479737999861181111857448896109510951094117495551210591079967588984747547796894584911101174979858795957661197710449879610119525109511109597944812677457117769104578968875710755511789967751054984117651049576691199910841071112997917844698612101112555499791091077664115975119796104117107747484104996899971297597104774669948659965411910978911747777996659566274675479458897105798947129595791051094688497510711591079796694771448911105984981276657910727999881046546117576756107699676747861047114476864611681010557957551055491057122969611861087969958549999128958927117911996965441010979976564555124671189546581079691068911464665977691075697699947105108895911661277867101069589712599551059712744679599467956978411868912455579612710694795107498879294969107991077486965296551158897481111811510976510797858101087461111125895711710791189456119397759999454481195995678565710596118991177610894106101111969779996968810116665410995647567558451099712997115977746441099577741010977575135610855555449596678769748811541057941010498105665981274115579755991155794611499115487945547510912459763699599569499667569878645105641097991051098899941085101210109116496951051179591211510751049656129746117710510671148591110105131046479277981149778745910978844791010910654511551058710879111097594109465474475549610101176575121010610411645105997561196510869499946394955102981067512118796558991149998756685612591154656104981099699767610677749497251157991011965571176910958485109791164845444861091054649910489576511957765955764478951164559897995455876116118989659541064691196997118976547104658776951097456644910557775794106410446576910105596951112698127891211125611557710677977104776966589445557524961161067896986749995994685585955119125988796884994871111106491010936114698959789961049584116911995576661112610777119479712811777559377410765712910968499785579107111012647644491251212510511868948567119976961286511986991011694710811109991011779249856810781066994677971010994477999669811994412481281041279765546467475696975767799578999987851061081056795910761046912696579697974976355105129898589779792951159911119967118549585472469596811810976651244641177599969995510811811410587979799109699494999787968105991098696898979999866979584947876710865647991269115954794696118588859842965485778994107549896846610544895797965112499107778487979116977766976449748461297499764955551244454955957451184589687510109879912799127957101299547910575107974965921065699874956576999469997549711899899667978488774109425117751146981995448956910689681295655951098710924776844479961065796544771081151154469795869107997469841161175711119796787991011961097127109468799551051061076795457862410971199799124710811987812912597889996945591059759511810755699896109101011698105774597663865955895744446710910859486666424945661069777467954510789847892485977785795969857104510791286422396789969757971077579377711547101069791048659410610974471158561081211048695789954798849964594749105747957777928967711787910127595610

Cellular component

cc_res <- res[rep(1:nrow(res), sapply(cc_class, length)), ]
cc_res$class <- unlist(cc_class)
library(metacoder)
data <- parse_taxonomy_table(input = as.data.frame(cc_res),
                             taxon_col = c("class" = -1),
                             other_col_type = "obs_info",
                             class_sep = "\\|")
data$obs_data$padj <- as.numeric(data$obs_data$padj)
data$taxon_funcs <- c(data$taxon_funcs,
                      change = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 obs_change <- data$obs_data[i, ]$log2FoldChange[data$obs_data[i, ]$padj <= min_p_value]
                                 mean(obs_change, na.rm = TRUE)
                               },
                               numeric(1))
                      },
                      num_changed = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 sum(data$obs_data[i, ]$padj <= min_p_value, na.rm = TRUE)
                               },
                               numeric(1))
                      })
set.seed(3)
data %>%
  filter_taxa(num_changed > 0) %>%
  filter_taxa(n_supertaxa <= 4) %>%
  # filter_taxa(n_supertaxa >= 1) %>% 
  filter_taxa(nchar(name) <= 40) %>%
  heat_tree(node_label = ifelse(nchar(name) <= 40, name, NA),
            node_size = num_changed,
            # node_size_trans = "log10",
            node_size_range = c(0.01, 0.03),
            # node_label_size_trans = "log10",
            node_label_size_range = c(0.01, 0.015),
            # edge_size_trans = "log10",
            edge_size_range = c(0.008, 0.03) / 2,
            node_color = 2^abs(change) * sign(change),
            node_color_trans = "linear",
            node_color_range = diverging_palette(),
            node_color_interval = c(-4, 4),
            edge_color_trans = "linear",
            edge_color_range = diverging_palette(),
            edge_color_interval =  c(-4, 4),
            node_label_max = 500,
            node_color_axis_label = "Factor change",
            node_size_axis_label = "Number of genes",
            layout = "da", initial_layout = "re",
            output_file = result_path("gene_expression--cellular_component"))

Biological Process

bp_res <- res[rep(1:nrow(res), sapply(bp_class, length)), ]
bp_res$class <- unlist(bp_class)
library(metacoder)
data <- parse_taxonomy_table(input = as.data.frame(bp_res),
                             taxon_col = c("class" = -1),
                             other_col_type = "obs_info",
                             class_sep = "\\|")
data$obs_data$padj <- as.numeric(data$obs_data$padj)
data$taxon_funcs <- c(data$taxon_funcs,
                      change = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 obs_change <- data$obs_data[i, ]$log2FoldChange[data$obs_data[i, ]$padj <= min_p_value]
                                 mean(obs_change, na.rm = TRUE)
                               },
                               numeric(1))
                      },
                      num_changed = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 sum(data$obs_data[i, ]$padj <= min_p_value, na.rm = TRUE)
                               },
                               numeric(1))
                      })
set.seed(7) #2, 4, 5, 7*, 19
mgsub <- function(pattern, replacement, x, ...) { # from: http://stackoverflow.com/questions/15253954/replace-multiple-arguments-with-gsub
  if (length(pattern)!=length(replacement)) {
    stop("pattern and replacement do not have the same length.")
  }
  result <- x
  for (i in 1:length(pattern)) {
    result <- gsub(pattern[i], replacement[i], result, ...)
  }
  result
}

to_replace <- matrix(ncol = 2, byrow = TRUE,
                     c("regulation of growth", "",
                       "activation of innate immune response", "",
                       "system development", "",
                       "regulation of response to stimulus", "",
                       "lipid metabolic process", "",
                       "selenium compound metabolic process", "selenium metabolic process"
                       ))
output_path <- file.path(output_folder,
                                    paste0("gene_expression--biological_process",
                                           output_format))
data %>%
  filter_taxa(num_changed > 0) %>%
  filter_taxa(n_supertaxa <= 3) %>%
  mutate_taxa(name = gsub("_", " ", name),
              f_change = 2^abs(change) * sign(change)) %>%
  mutate_taxa(short_name = vapply(name, FUN.VALUE = character(1), function(x) {
    mgsub(pattern = to_replace[, 1], replacement =  to_replace[, 2], x, fixed = TRUE)
  })) %>%
  heat_tree(node_label = ifelse(abs(f_change) > 1, short_name, NA),
            node_size = num_changed,
            # node_size_trans = "log10",
            node_size_range = c(0.008, 0.03),
            # node_label_size_trans = "log10",
            node_label_size_range = c(0.012, 0.02),
            # edge_size_trans = "log10",
            edge_size_range = c(0.008, 0.03) / 2,
            node_color = f_change,
            node_color_trans = "linear",
            node_color_range = diverging_palette(),
            node_color_interval = c(-5, 5),
            edge_color_trans = "linear",
            edge_color_range = diverging_palette(),
            edge_color_interval =  c(-5, 5),
            node_label_max = 80,
            node_color_axis_label = "Fold change",
            node_size_axis_label = "Number of genes",
            layout = "da", initial_layout = "re",
            output_file = result_path("gene_expression--biological_process"))

Molecular Function

mf_res <- res[rep(1:nrow(res), sapply(mf_class, length)), ]
mf_res$class <- unlist(mf_class)
library(metacoder)
data <- parse_taxonomy_table(input = as.data.frame(mf_res),
                             taxon_col = c("class" = -1),
                             other_col_type = "obs_info",
                             class_sep = "\\|")
data$obs_data$padj <- as.numeric(data$obs_data$padj)
data$taxon_funcs <- c(data$taxon_funcs,
                      change = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 obs_change <- data$obs_data[i, ]$log2FoldChange[data$obs_data[i, ]$padj <= min_p_value]
                                 mean(obs_change, na.rm = TRUE)
                               },
                               numeric(1))
                      },
                      num_changed = function(x, subset = NULL) {
                        vapply(obs(x),
                               function(i) {
                                 sum(data$obs_data[i, ]$padj <= min_p_value, na.rm = TRUE)
                               },
                               numeric(1))
                      })
data %>%
  filter_taxa(num_changed > 0) %>%
  filter_taxa(n_supertaxa <= 3) %>%
  # filter_taxa(n_supertaxa >= 1) %>% 
  filter_taxa(nchar(name) <= 40) %>%
  heat_tree(node_label = ifelse(nchar(name) <= 40, name, NA),
            node_size = num_changed,
            # node_size_trans = "log10",
            node_size_range = c(0.01, 0.03),
            # node_label_size_trans = "log10",
            node_label_size_range = c(0.01, 0.015),
            # edge_size_trans = "log10",
            edge_size_range = c(0.008, 0.03) / 2,
            node_color = 2^abs(change) * sign(change),
            node_color_trans = "linear",
            node_color_range = diverging_palette(),
            node_color_interval = c(-4, 4),
            edge_color_trans = "linear",
            edge_color_range = diverging_palette(),
            edge_color_interval =  c(-4, 4),
            node_label_max = 500,
            node_color_axis_label = "Factor change",
            node_size_axis_label = "Number of genes",
            layout = "da", initial_layout = "re",
            output_file = result_path("gene_expression--molecular_function"))

Software and packages used

sessionInfo()
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.2 LTS
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] stats4    parallel  grid      stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] metacoder_0.1.2            org.Hs.eg.db_3.3.0        
##  [3] GO.db_3.3.0                AnnotationDbi_1.34.4      
##  [5] DESeq2_1.12.4              airway_0.106.2            
##  [7] SummarizedExperiment_1.2.3 Biobase_2.32.0            
##  [9] GenomicRanges_1.24.3       GenomeInfoDb_1.8.7        
## [11] IRanges_2.6.1              S4Vectors_0.10.3          
## [13] BiocGenerics_0.18.0        rnaseqGene_0.99.121218    
## [15] knitcitations_1.0.7        knitr_1.14                
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-128                  bitops_1.0-6                 
##  [3] matrixStats_0.51.0            lubridate_1.6.0              
##  [5] RColorBrewer_1.1-2            httr_1.2.1                   
##  [7] rprojroot_1.2                 fission_0.106.2              
##  [9] tools_3.3.1                   backports_1.0.5              
## [11] R6_2.2.0                      rpart_4.1-10                 
## [13] Hmisc_3.17-4                  DBI_0.5-1                    
## [15] lazyeval_0.2.0                Gviz_1.16.5                  
## [17] mgcv_1.8-15                   colorspace_1.2-7             
## [19] nnet_7.3-12                   gridExtra_2.2.1              
## [21] formatR_1.4                   labeling_0.3                 
## [23] rtracklayer_1.32.2            scales_0.4.1                 
## [25] genefilter_1.54.2             stringr_1.1.0                
## [27] digest_0.6.12                 Rsamtools_1.24.0             
## [29] foreign_0.8-67                rmarkdown_1.3                
## [31] XVector_0.12.1                dichromat_2.0-0              
## [33] htmltools_0.3.5               bibtex_0.4.0                 
## [35] PoiClaClu_1.0.2               ensembldb_1.4.7              
## [37] BSgenome_1.40.1               RSQLite_1.0.0                
## [39] BiocInstaller_1.22.3          shiny_0.14.1                 
## [41] BiocParallel_1.6.6            dplyr_0.5.0                  
## [43] acepack_1.3-3.3               VariantAnnotation_1.18.7     
## [45] RCurl_1.95-4.8                magrittr_1.5                 
## [47] Formula_1.2-1                 Matrix_1.2-7.1               
## [49] Rcpp_0.12.9                   munsell_0.4.3                
## [51] RefManageR_0.13.1             stringi_1.1.2                
## [53] yaml_2.1.13                   RJSONIO_1.3-0                
## [55] zlibbioc_1.18.0               plyr_1.8.4                   
## [57] AnnotationHub_2.4.2           lattice_0.20-34              
## [59] Biostrings_2.40.2             splines_3.3.1                
## [61] GenomicFeatures_1.24.5        annotate_1.50.1              
## [63] locfit_1.5-9.1                igraph_1.0.1                 
## [65] reshape2_1.4.2                geneplotter_1.50.0           
## [67] biomaRt_2.28.0                XML_3.98-1.4                 
## [69] evaluate_0.10                 biovizBase_1.20.0            
## [71] latticeExtra_0.6-28           data.table_1.10.4            
## [73] httpuv_1.3.3                  gtable_0.2.0                 
## [75] assertthat_0.1                ggplot2_2.2.1                
## [77] mime_0.5                      xtable_1.8-2                 
## [79] survival_2.39-5               tibble_1.2                   
## [81] pheatmap_1.0.8                GenomicAlignments_1.8.4      
## [83] cluster_2.0.4                 sva_3.20.0                   
## [85] interactiveDisplayBase_1.10.3 BiocStyle_2.0.3

Comments